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We present an efficient Monte-Carlo method for long-range interacting systems to calculate 
free energy as a function of an order parameter. In this method, a variant of the Wang-Landau 
method regarding the order parameter is combined with the stochastic cutoff method, which has 
recently been developed for long-range interacting systems. This method enables us to calculate 
free energy in long-range interacting systems with reasonable computational time despite the 
fact that no approximation is involved. This method is applied to a three-dimensional magnetic 
dipolar system to measure free energy as a function of magnetization. By using the present 
method, we can calculate free energy for a large system size of 16 3 spins despite the presence of 
long-range magnetic dipolar interactions. We also discuss the merits and demerits of the present 
method in comparison with the conventional Wang-Landau method in which free energy is 
calculated from the joint density of states of energy and order parameter. 

KEYWORDS: Monte Carlo, long-range interacting system, Wang-Landau method, free energy measure- 
ment, magnetic dipolar system 



In general, Monte Carlo (MC) simulations in long- 
range interacting systems are much more difficult than 
those in short-range interacting systems because we have 
to take a large number of interactions into consideration. 
For example, in the case of systems with pairwise in- 
teractions, the number of interactions is proportional to 
N 2 , where N is the number of elements of the system. 
Therefore, if a naive MC simulation is carried out in such 
systems, the computational time per MC step rapidly in- 
creases in proportion to N 2 , which is in contrast to that 
in the case of short-range interacting systems in which 
the computational time increases in proportion to N. In 
order to overcome this difficulty, many simulation meth- 
ods have been proposed. 1 ~ 11 ) 

Recently, one of the authors and Matsubara have de- 
veloped an efficient MC method called the stochastic cut- 
off (SCO) method for long-range interacting systems. 12 ) 
In the SCO method, each of the pairwise interactions 
Vij is stochastically switched to cither or a pscudoint- 
eraction Vy by the stochastic potential switching (SPS) 
algorithm. 13 ' 14 ^ The switching probability to and that 
to Vij are Pij and 1 — Py, respectively. Since the pseu- 
dointeraction Vij and switching probability Py are cho- 
sen properly in the SPS algorithm, the SCO method 
strictly satisfies the detailed balance condition concern- 
ing the original Hamiltonian. 13-15 ) This means that the 
SCO method does not involve any approximation. Fur- 
thermore, since most of the distant and weak interac- 
tions are switched to and an efficient method to switch 
potentials has been developed, 12 ) the SCO method en- 
ables us to reduce the computational time of long-range 
interactions markedly. For example, in the case of three- 
dimensional dipolar systems, to which our new method 
will be applied later, the computational time is reduced 
from 0(N 2 ) to 0(N log N) by the SCO method. 12 ) We 
can measure internal energy and heat capacity without 
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0(N ) computation by a method proposed in ref. 15. 
An efficient method to combine the SCO method with 
the replica exchange method 16 ) has also been developed 
there. 

In this letter, we propose an efficient MC method 
of combining the SCO method with the Wang-Landau 
method. 17,18 ) This method enables us to calculate free 
energy as a function of an order parameter with rea- 
sonable computational time even in long-range interact- 
ing systems. In the case of three-dimensional magnetic 
dipolar systems, computational time per MC step is re- 
duced from 0{N 2 ) to 0(N log N). As will be shown 
later, this method enables us to calculate free energy in 
a three-dimensional magnetic dipolar system with a size 
of 16 3 spins. This system size is much larger than that of 
the previous work, 19 ) i.e., 10 3 spins, in spite of the fact 
that long-range dipolar interactions are included in the 
present work and they are not included in the previous 
work. 

We now start to present our new MC method. As an 
example, we hereafter consider to measure free energy as 
a function of the z-component of the magnetization in 
a classical Hcisenberg spin system. It is straightforward 
to generalize the method for other cases. Free energy is 
defined by 



cxp[-/3P(/3;m z )] 

= CTt {Sf} exp[-pH{Si}]6(r 



where 



m*{Si} 



l{Si}), (1) 



(2) 



where TV is the number of spins, j3 is the inverse tem- 
perature, T~L{Si} is the Hamiltonian of the spin system, 
and C is a constant. It is not important how we choose 
C because it only contributes to F((3; m z ) as a constant. 
The right-hand side of eq. (1) is the sum of the weights 
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of all the states with a magnetization m z . 

In the present method, we use a variant of the Wang- 
Landau method regarding the order parameter. A similar 
method has been proposed in ref. 20 to calculate free en- 
ergy by the multicanonical ensemble method. 21,22 ' The 
basic idea is as follows: We perform an MC simulation 
with the Hamiltonian 

H'{S i } = H{S i }-p- 1 G(m* z {S i }). (3) 

During the simulation, G(m z ) in eq. (3) is modified so 
that P(m z ) becomes a constant, where P{m z ) is the 
probability that a state with a magnetization m z is sam- 
pled. Then, the resultant function G(m z ) is related to 
the free energy F([3;m z ) by 

G(m z ) = /3F(f3; m z ) + constant. (4) 

This can be easily shown as 

P(m z ) 

oc Tr^SV} exp[-/3'H'{Si}]S(m z - m* z {Si}) 

= exp[G(m,)]Tr { £ !} cxp{- f3n{S t }]d(m z - m* z {S t }) 

— cxp[G(m z )]C^ 1 cxp[— /3F(/3; m z )} = constant, (5) 

where we have used eq. (1) to go from the third line to 
the fourth. To modify G(m z ) so that P(m z ) becomes a 
constant, we use the conventional procedure of the Wang- 
Landau method, 17 ' 18 ^ i.e., we modify G(m z ) as 

G(m z ) G(m z ) - AF, (AF > 0) (6) 

after each trial to update a single spin. If we start our 
simulation with G(m z ) = 0, as conventionally performed, 
states with low free energies are frequently sampled 
at the beginning of the simulation. However, since the 
weights of such states are reduced more by the reduction 
in G{m z ) [recall that the Hamiltonian is given by eq. (3)], 
G(m z ) is adjusted by this procedure so that P(m z ) be- 
comes a constant. As is conventionally performed in the 
Wang-Landau method, 17, 18 ' the constant AF in eq. (6) 
is gradually reduced as the simulation proceeds. 

Now, we briefly compare the present method with the 
conventional Wang-Landau method in which one evalu- 
ates the joint density of states as a function of energy 
and magnetization 

n(E, m z ) ee Tt {Sz} S(E - H{Si})5(m z - m* z {S t }), (7) 
and calculates free energy from n(E, m z ) as 

/oo 
dEexp(-/3E)n(E,m z ). (8) 

In the conventional Wang-Landau method, we modify a 
two- variable function G(E,m z ), which becomes propor- 
tional to n(E,m z ) at the end of the simulation, by a 
method similar to eq. (6). Therefore, we have to calcu- 
late E and m z after each update of a single spin. Since 
the number of interactions per spin is N — 1 in long-range 
interacting systems, the computational time for calculat- 
ing the new energy is O(N). In contrast, we only need to 
calculate m z in the present method, and the computa- 
tional time to calculate the new magnetization is 0(1). 
This is the reason why we use the variant of the Wang- 



Landau method. 

Note that the above-mentioned simulation method is 
still time-consuming if % in eq. (3) involves long-range 
interactions because we need to carry out an MC simu- 
lation with the Hamiltonian H! . To overcome such dif- 
ficulty, we simply use the SCO method. Since the SCO 
method with the Hamiltonian TL' samples a state accord- 
ing to the Boltzmann weight cxp(— fiH.') as in the conven- 
tional MC method, we can measure F(/3; m z ) in the same 
way as before. In the case of three-dimensional magnetic 
dipolar systems, the number of interactions per spin is 
reduced from N - 1 to £>(log N) by the SCO method. 12 ' 
Therefore, the computational time per MC step becomes 
0(N log N). Because the SCO method can be applied to 
a part of the Hamiltonian, 13,15 ' we only apply the SCO 
method to a long-range part in %. 

To summarize, we show the whole procedure of our 
method. When the original Hamiltonian consists of long- 
range interactions H^{Si} = J2i<j Sj) {Vij i s a 
pairwise interaction) and short-range ones w^' '{Si}, the 
procedure proceeds as follows: 

1) Set G(m z ) = and H{m z ) = 0, where H(m z ) is a 
histogram to check whether or not all the magne- 
tizations are sampled with equal probabilities. The 
initial AF in eq. (6) is set sufficiently large so that 
G(m z ) is adjusted quickly in the early stage of the 
Wang-Landau method. 

2) Repeat the following two steps as a basic MC pro- 
cedure: 

a) Switch each of Vij in %( L ' to cither or Vj with 
a probability of P^ or 1 — P^j, respectively. The 
method proposed in ref. 12 is used to switch po- 
tentials efficiently. The probability Pij is 

PijiSt, Sj) = cxp[/3(^-(S 4 , Sj) - V*)], (9) 

where /? is the inverse temperature and V*j is a 
constant equal to (or greater than) the maximum 
Vij (Si, Sj). The pseudopotential Vij is defined by 

VijiS^Sj) ee VijiSi, S^-r 1 log[l-P y (S i; Sj)}. 

(10) 

This potential switching step is performed every 
Switch MC steps. 

b) Perform a standard MC simulation with the 
Hamiltonian 

H'{Si} = H^iSij+Y^VijiS^SjyGimliSi}), 

(11) 

where the sum ^ runs over potentials which are 
switched to Vij in step a). During the simulation, 
G(m z ) is changed according to eq. (6) after each 
trial to update a spin. This adjustment of G(m z ) 
is performed regardless of whether or not the trial 
is accepted. We also change the histogram as 

H{m z ) ->> H(m z ) + 1, (12) 

after each trial. 

3) Check whether or not the histogram H(m z ) is flat. 
If it is flat, halve AF and reinitialize the histogram 
as H(m z ) = 0. This check is performed every £ c hcck 
MC steps. 



J. Phys. Soc. Jpn. Letter. Online- Journal Subcommittee 3 




Fig. 1. (Color online) m z dependences of F(/3;m z ) for D = 0, 

0.02, 0.04, 0.06, 0.08, and 0.1 J. The system size N is 16 3 and Fig- 2. (Color online) Average number n of potentials per site 

the temperature T is 0.7 J. All the data are shown by lines that arc switched to Vij ^ is plotted as a function of D for 

without error bars. Symbols are drawn with error bars at several T = 0.7, 0.9, 1.1, and 1.3 J (from top to bottom), 
data points. The average is taken over 10 different runs. 



4) Stop the simulation if AF is small enough. Other- 
wise, return to 2). 

To check the efficiency of the method, we apply it 
to measure free energy as a function of m z in a three- 
dimensional magnetic dipolar system. The Hamiltonian 
is given as 



-jy St 



z -.2 



Sj Cu / 



Sj ■ Sj ^ (Sj ■ Tij)(Sj ■ Tjj) 



,(13) 



where Si is a classical Hciscnbcrg spin of \Si\ = 1, 
(ij) runs over all the nearest-neighboring pairs, is 
a vector spanned from site i to site j in the unit of the 
lattice constant a, and ry = \rij\- On the right-hand 
side of eq. (13), the first term describes ferromagnetic 
exchange interactions, the second term uniaxial magne- 
tocrystalline anisotropy energies whose easy axis is paral- 
lel to the z-direction, and the third term magnetic dipo- 
lar interactions. As mentioned above, the system size N 
is 16 3 . The boundary condition is open in all directions. 
We fix the ratio C u / J to 0.1 and change D/ J to see how 
the structure of free energy depends on the strength of 
dipolar interactions. 

In the present simulations, we set the initial AF in 
eq. (6) to J. We stop our simulation after we halve AF 
20 times. Therefore, the final AF is J x 2~ 20 . We have 
checked that G(m z ) converges well in later stages of the 
Wang-Landau method. The histogram H(m z ) is checked 
every 10,000 MC steps. We regard the histogram as flat 
when H(m z ) for all the magnetizations is not less than 
80% of the average histogram (H(m z )). We estimate 
F(/3;m z ) in the range of < m z < 0.98 on a grid of 
40,142 bins. Since F(/3;m z ) is an even function of m z , 
we calculate F(/3;m z ) for positive magnetizations. The 
SCO method is applied only to dipolar interactions. They 
are switched every 10 MC steps. V*j in eq. (9) is set to 
the maximum Vij (Si, Sj). 



The result of the free energy measurement at T = 0.7 J 
is shown in Fig. 1, where we set the Boltzmann constant 
fee to unity. The temperature is well below the critical 
temperature of the model for C u = D = 0, which is es- 
timated to be about 1.44 J. 23 ) For each D, we carried 
out 10 different runs with different initial conditions and 
random sequences to estimate the means and error bars 
of the data. As a result, we have found that the error 
bars are less than 1.0 J for all the data. In Fig. 1, sym- 
bols are drawn with error bars. The error bars are much 
smaller than the symbols. From the smallncss of the er- 
ror bars, we consider that correct data are obtained by 
the present method. We also see that the position of the 
global minimum changes from m z ks 0.8 to m z = as D 
increases. This result is reasonable because dipolar inter- 
actions prefer demagnetized states. The case D = 0.04 J 
is a marginal one in which free energies at the two min- 
ima are almost the same. 

To estimate the efficiency of the present method, we 
first measure the average number n of potentials per site 
that survive as Vij for several temperatures and D's. The 
result is shown in Fig. 2. The average number n increases 
with decreasing temperature and increasing D. However, 
it is about 7 even when T = 0.7 J and D — 0.1 J. 
This means that more than 99.8% of the interactions are 
cut off by being switched to Vij = 0. We next examine 
how £mc depends on temperature, where £mc is the to- 
tal number of MC steps until AF is halved 20 times and 
the simulation is stopped. Figure 3 shows the result. The 
average is taken over 50 different runs. We measure £mc 
for D/J = and D/J = 0.04. We do not use the SCO 
method in the former case because long-range dipolar in- 
teractions are absent. In both cases, £mc increases with 
decreasing temperature since relaxation becomes slower 
at lower temperatures. £mc f° r D/J~ 0.04 increases 
more rapidly than that for D/J = 0. However, the tem- 
perature dependence is not so strong. We also find that 
^mc's for D/J = and D/J = 0.04 are not very differ- 
ent. This result shows that £mc does not increase much 
by the use of the SCO method. The computational time 
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Fig. 3. (Color online) Temperature dependences of the total num- 
ber of MC steps £mc m the free energy calculation for D = 
(full circle) and D = 0.04 J (full square). The average is taken 
over 50 different runs. 



per run for T = 0.7 J and D = 0.1 J, which is the most 
time-consuming case we have examined, was about ten 
days by a single-core calculation with a Core-i7 2.8 GHz 
processor. 

Finally, we again compare the present method with the 
conventional Wang-Landau method in which the joint 
density of states defined by cq. (7) is evaluated. The 
first merit of the present method is that the computa- 
tional time for long-range interactions is markedly re- 
duced by the use of the variant of the Wang-Landau 
method and the SCO method. As mentioned before, the 
difficulty in the conventional Wang-Landau method in 
long-range interacting systems is that we have to calcu- 
late the Hamiltonian H{Si} after each update of a single 
spin. It is very desirable to develop a method of com- 
bining the conventional Wang-Landau method with the 
SCO method. To this end, an approach used in ref. 15 
might be helpful. The second merit of the present method 
is that the function G(m z ) to be adjusted is a one vari- 
able function. In contrast, we have to adjust a two- 
variable function G(E,m z ) in the conventional Wang- 
Landau method. This is the main reason why the system 
size accessible by the present method is larger than that 
by the conventional Wang-Landau method. 19 ^ However, 
the trade-off for this merit is that the temperature is 
kept constant. Therefore, in the present method, we can 
only estimate free energy at one temperature by a single 
simulation. In contrast, we can estimate free energy at 
any temperatures by a single simulation of the conven- 
tional Wang-Landau method because free energy at any 
temperatures can be calculated from the joint density 
of states using eq. (8). Furthermore, the present method 
has another drawback when the simulation is performed 
at low temperatures. In the conventional Wang-Landau 
method, high-energy states with high entropies are the 
source of fast relaxation, and the system rapidly forgets 
the present state when the system reaches a high-energy 
region. However, no such source exists in the present 



method when the temperature is low. Note that, in the 
present method, the zero-magnetization state is not a 
high-entropy state. When the temperature is low, the 
present method samples only a small portion of the states 
with low energies at any magnetization. This means that 
it is not trivial in the present method that equilibrium 
sampling is realized at low temperatures. Therefore, we 
should carefully check whether or not equilibrium sam- 
pling is realized. As performed in the present work, an 
effective way to check equilibration is by measuring free 
energy several times using different initial conditions and 
random sequences and by checking whether or not the 
same result is obtained. One should keep in mind that 
the present method has these drawbacks. 

In summary, we have developed an efficient MC 
method of free energy calculation in long-range inter- 
acting systems by combining a variant of the Wang- 
Landau method with the stochastic cutoff method. The 
efficiency of the method has been confirmed by ap- 
plying the method to a free energy calculation in a 
three-dimensional magnetic dipolar system. We have 
also discussed the merits and demerits of the present 
method in comparison with the conventional Wang- 
Landau method. 

The authors would like to thank Professor K. Sasaki 
for valuable discussions and comments. This work is sup- 
ported by a Grant-in-Aid for Scientific Research (No. 
21740279) from MEXT. 
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